Tension zone trapped by exogenous cline: Analysis of a narrow hybrid zone between two parapatric Oxytropis species (Fabaceae)

Abstract Hybrid zones have been widely highlighted for their interest in understanding evolutionary processes. It is generally accepted that hybrid zones can be maintained in a balance between dispersal and selection. However, the selective forces can either be endogenous (i.e., genetic incompatibilities between parental taxa) or exogenous (i.e., parental taxa are adapted to different environments). In this study, we evaluated these alternatives and determined the maintenance of a narrow hybrid zone between parapatric distributed Oxytropis diversifolia and O. leptophylla in Nei Mongol, China. For 507 individuals sampled from two populations in the hybrid zone, 12 O. diversifolia populations and five O. leptophylla populations, we measured leaf‐morphological characteristics, quantified genetic structure using 11 nuclear microsatellite loci and five chloroplast DNA intergenic regions, collected micro‐ and macrohabitat data, and conducted geographical cline analysis. We found that the two species differed in leaf morphology, and putative hybrids showed either intermediacy or a bias to O. diversifolia. Parental taxa formed two genetically distinct clusters, while populations in the hybrid zone consisted of both parental forms and various admixed individuals, exhibiting a bimodal pattern. The hybrid zone was coupled to ecological transitions of both microhabitat (i.e., the slope) and macroclimatic conditions. However, the genetic clines were significantly narrower than the environmental cline. Our results indicate that endogenous selection can be primarily responsible for maintaining the hybrid zone, while local adaptation accounts for the position of the zone. We further suggest the probable outcome of hybridization could be introgression.


| INTRODUC TI ON
Natural hybridization, the mating of individuals from distinguishable populations, usually of distinct taxa, frequently occurs in plants and can have diverse evolutionary consequences (Abbott et al., 2013;Buerkle et al., 2003;Harrison, 1993;Hopper, 1995;Whitney et al., 2010). Hybridization may slow or reverse differentiation by allowing gene flow and recombination, and consequently results in breakdown of species boundaries and cause extinction of rare endemics (Buerkle et al., 2003;Levin et al., 1996;Rhymer & Simberloff, 1996;Wolf et al., 2001). On the other hand, hybridization may contribute to speciation through the formation of new hybrid taxa (either via allopolyploidy or homoploid hybrid speciation), and introgression of a few loci may also promote adaptive divergence and so facilitate speciation (reviewed in Abbott et al., 2013).
Moreover, divergent lineages with incomplete barriers of reproduction may form stable hybrid zones, which can be examined to determine mechanisms that reduce gene flow and genome merging (Abbott, 2017;Stankowski et al., 2021). Depending on the hybridizing taxa involved, an assessment of the evolutionary outcome of hybridization can provide important insights in predicting species persistence and hybrid speciation, as well as the maintenance of biological diversity (Yan et al., 2017).
Most hybrid zones have existed for an extended time, so they integrate the outcome of evolutionary processes over timescales that are not accessible in laboratory experiments (Nolte et al., 2006).
Hybrid zones often originate from divergent lineages that are evolving independently, and come into secondary contact at the edges of their distributions. The classic hybrid zone theory (Fisher, 1950;Haldane, 1948) predicts that the extent and shape of hybrid zones depend upon a balance between migration and selection, and the latter often involves a combination of "endogenous" and "exogenous" forms of selection (Barton, 2001;Moore & Price, 1993). Endogenous selection expects that internal genetic incompatibilities reduce the viability or reproductive success of hybrids, and hybrid zones are maintained by a balance between random dispersal and selection against hybrids, which is so-called "tension zone model" (Barton & Hewitt, 1985). In this model, the antihybridization mechanism would be expected independent of ecological conditions (Swenson, 2008).
Exogenous selection assumes that the fitness of genotypes is coupled to the environment, and that selection favors different genotypes in different environments (Barton & Gale, 1993;Moore & Price, 1993). As a result, hybrid zones tend to be located in ecotones between habitats occupied by parental forms (i.e., ecotone model; Nolte et al., 2006). An ecotone model predicts selection against pure genotypes in alternative habitats, and less selection against hybrids if the latter have an intermediate degree of habitat dependent fitness (Nolte et al., 2006).
The hybrid zone literature has long debated the relative roles of endogenous and exogenous selection (Arnold, 1997;Barton, 2001), and in many cases both exogenous and endogenous factors are acting concomitantly (Bierne et al., 2011). One way to evaluate these alternatives is to do geographical cline analysis, comparing the width of genetic cline with that of the environmental cline (e.g., Stankowski, 2013;Yanchukov et al., 2006). Furthermore, theoretical modeling and literature review suggest that tension zones have a tendency to become trapped by, and therefore coincide with, exogenous barriers due to ecological selection (Bierne et al., 2011).
The result is that local adaptation explains where genetic breaks are positioned, but not necessarily their existence, which can be best explained by endogenous incompatibilities (Bierne et al., 2011).
Additional insight can be gained by taking a multilocus approach and examining the frequency distribution of genotypic classes within a hybrid zone. Where a hybrid zone is maintained by weak to moderate postzygotic selection rather than prezygotic barriers, a unimodal distribution is predicted; where postzygotic selection is strong, or prezygotic barriers are present, the distribution is expected to become flatter and, ultimately, bimodal (Jiggins & Mallet, 2000).
It is generally accepted that closely related species tend to hybridize more often. In particular, species in rapidly diversifying adaptive radiations may be particularly prone to hybridization (Abbott et al., 2013;Gourbière & Mallet, 2010;Seehausen, 2004).
The genus Oxytropis DC. (Fabales: Fabaceae) is recognized to have a history of recent rapid radiation (Kholina et al., 2016;Shahi Shavvon et al., 2017). Oxytropis is comprised of about 310 species mainly distributed in temperate regions of the Northern Hemisphere (Zhu et al., 2010), among which the vast majority are restricted to Central and East Asia (Zhao & Liu, 1997). Using nuclear (ITS) and plastid (trnL-F) sequence data, Shahi Shavvon et al. (2017) found that Oxytropis splits from the close relative Astragalus at ca. 15.6 million years ago in the mid Miocene, and the crown group of Oxytropis was estimated to date back to 5.6 Ma in the late Miocene. A recent study, which employed the chloroplast genome sequences based on the genome skimming approach, revealed that Oxytropis is sister to the Coluteoid clade (Astragalus is sister to the Oxytropis + Coluteoid clade), and the divergence of Oxytropis from the Coluteoid clade was dated to 14.08 Ma (Su et al., 2021). It is suggested that hybridization might have played a significant role in the evolution of Oxytropis F I G U R E 1 (a) Map of northern China showing locations of 19 sampled Oxytropis populations. Each pie chart indicates the proportions of different leaf-shape phenotypes in the population. Twelve populations are O. diversifolia; Δ two populations in the hybrid zone (G and HUA); *five populations of O. leptophylla. The dashed line is the contour where the alternative cpDNA haplotypes are predicted to be at a frequency of 0.5. The small inset shows occurrence records of two parental species and the hybrid zone. See section 2.7.2 for details. White circles, O. diversifolia; black circles, O. leptophylla; gray circles, two populations in the hybrid zone. (b-g) phenotypes of Oxytropis diversifolia (a, 1 leaflet; b, 1-3 leaflets; c, 3 leaflets. Flower color: Pale yellow), O. leptophylla (d, 5-13 leaflets. Flower color: Purplish red) and putative hybrids (e, 3-5 leaflets with pale yellow flowers; f, 1-3 leaflets with intermediate flower color). Scale bars = 5 mm. (Kholina et al., 2016;Malyshev, 2008;Shahi Shavvon et al., 2017).
However, the evidence of natural hybridization is comparatively rare, and is based on either morphology or plastid sequence data (e.g., Kholina et al., 2019Kholina et al., , 2021. The consequence of natural hybridization has not been explored yet. In this study, we analyze a secondary contact between parapatric Oxytropis diversifolia E. Peter (Fabales: Fabaceae) and O. leptophylla (Pall.) DC. in Nei Mongol, China. Both species belong to sect.
Xerobia, which is comprised of 27 species (Zhu et al., 2010) that are mostly cryoxerophytes associated with mountain-steppe territories (Kholina et al., 2021). O. diversifolia is mainly distributed in desertsteppe regions of Nei Mongol, whereas O. leptophylla can be found in typical-steppe and forest-steppe regions, and is more widely distributed in northeastern China (Nei Mongol, Jilin, Hebei, Shanxi), Mongolia and Russia (Zhao & Liu, 1997;Zhu et al., 2010;Figure 1a). The clinal pattern of the heritable leaf-shape variation can be attribute to natural selection imposed by spatially varying macroclimatic factors, other than to neutral genetic structure (Wang et al., 2021). (Figure 1e; Table 1), and the main differences include imparipinnate leaves with 5-13 leaflets and a purple corolla.
The distributions of these two species overlap close to the longitude of 110° E, where we found two populations ( Figure 1a) in which individuals with intermediate morphology occur (e.g., Figure 1f,g; for details of the morphological differences, see Table 1). It promoted us using molecular markers to verify whether these peculiar individuals represent the offspring of natural hybridization between these species, and if so, to explore what the consequence of natural hybridization is. Furthermore, the hybrid zone is located in ecotones between desert-steppe regions of O. diversifolia and typical-steppe regions of O. leptophylla, and therefore we would like to test whether exogenous selection (i.e., microhabitat and/or macrohabitat factors) against dispersing individuals contribute to shaping and maintaining the hybrid zone.
Here, we use an integrative approach by combining morphological measurements, molecular data, and environmental association analyses to evaluate the pattern and maintenance of the hybrid zone between O. diversifolia and O. leptophylla. Specifically, we seek to address the following questions: (1) Considering leaf-morphological traits, what is the phenotypic difference between parental species and putative hybrids? (2) What is the genetic structure of the hy- Zhao (Wang et al., 2021). O. leptophylla and O. diversifolia have the same number of chromosomes (2n = 16) (Zhang et al., 1994). O. diversifolia is flowering from April to May, while the flowering of O. leptophylla occurs slightly later from May to June (Zhu et al., 2010; H. Wang, personal observation). Because the two species and putative hybrids in the hybrid zone have quite similar floral structure (papilionaceous corolla) and floral size (Zhu et al., 2010; personal observation), we mainly measured leaf-morphological traits in this study. Purplish red, rarely blue a Putative hybrids from two populations in the hybrid zone were identified based on nuclear microsatellite data, and may include F 1 /F 2 hybrids, backcross types and various segregation products.

| Population survey and sampling
A total of 19 populations were surveyed in April and May from 2016 to 2018 ( Figure 1a; Table S1). Based on leaf-shape phenotype (i.e., the number of leaflets) and flower color, localities correspond to putatively "pure" populations of O. diversifolia (N = 12, covering its entire distribution range in Nei Mongol), O. leptophylla (N = 5, representatives in eastern Nei Mongol), or of the hybrid zone with signs of intermediacy (N = 2). In each population, leaf-phenotype survey was performed by tallying leaf-phenotype counts (≤100 individuals), or along three to five transects (3 m wide) spanning the length of the population (>100 individuals). This survey was used to estimate population size (i.e., the total number of individuals) and population density (population size divided by total area). Silica-gel dried leaf samples were collected, using methods described in Wang et al. (2021). In total, 507 individuals were col-

| Nuclear microsatellite analyses
2.5.1 | Population genetics GENEPOP 4.7 (Rousset, 2008) was used to test for genotypic linkage disequilibrium between each pair of loci within each population, and the probability was adjusted by Benjamini-Hochberg correction (Benjamini & Hochberg, 1995). An assessment of departure from Hardy-Weinberg equilibrium was conducted using exact tests (Rousset & Raymond, 1995) implemented by GENEPOP. For each population, mean allele number per locus (N a ), observed heterozygosity (H o ), unbiased expected heterozygosity (H e ) and F IS , were calculated using GenAlEx version 6.503 (Peakall & Smouse, 2006. Those genetic diversity parameters were compared to explore whether they differ between populations located outside and within the hybrid zone.

| Parent and hybrid assignment
To estimate the extent of admixture between parental species in the hybrid zone, the Bayesian clustering program STRUCTURE 2.3.4 (Pritchard et al., 2000) was used to assign all 507 individuals. An admixture model with independent allele frequencies were applied, and analyses were conducted with a K value that varied from 1 to 10. For each K value, analysis was repeated 10 times with a burn-in of 10,000 Markov Chain Monte Carlo (MCMC) steps followed by 50,000 iterations. The best-fit number of distinct groups (K) was determined based on the maximum ΔK value (Evanno et al., 2005) using CLUMPAK (Kopelman et al., 2015). Tq is a threshold value used to determine which class an individual should belong to (De hert et al., 2012). Tq values previously used to assign parental and hybrid individuals have varied from 0.80 to 0.95 (Yan et al., 2017).
In our STRUCTURE analysis, a relatively strict threshold Tq value of 0.90 was used, above which the majority of individuals from parental populations were correctly assigned ( Table S4). The Q scores of individuals in the hybrid zone were then binned in 0.10 increments and compiled in a histogram. The shape of the frequency distribution (either unimodal or non-unimodal) was determined by Hartigan's dip test (Hartigan, 1985;Hartigan & Hartigan, 1985).
To check the assignment of all individuals to parents (P 1 , P 2 ) and leptophylla. Second, NEWHYBRIDS Version 1.1 beta (Anderson & Thompson, 2002) was used to compute the posterior probability that individuals fall into six different purebred and hybrid categories.
The prior information of the parental populations was given (z0 = O. diversifolia and z1 = O. leptophylla). Additional parameter settings included a burn-in of 10,000 MCMC steps followed by 100,000 iterations. A relatively strict posterior probability of 0.95 was used.
Results of those two analyses were quite similar to STRUCTURE results, but discordance still existed ( Table S4). The final assignments were retrieved with the criteria that two out of the above three methods yielded consistent results.  (Nei, 1987). A statistical parsimony network based on Median Joining method (Bandelt et al., 1999) was constructed and visualized using PopArt 1.7 (Leigh & Bryant, 2015).

| Microhabitat differentiation
The microhabitat axes captured habitat attributes at the local scale of individuals within populations. For each individual sampled, a 40 cm × 40 cm square plot was surveyed ( Figure S1), slope (by geologic compass), the percent of total vegetation cover, the percent of rocky ground (rock diameter >0.2 cm), and the percent of bare ground (sand and soil) were measured as described in Wang et al. (2021). In total, microhabitat data for 477 individuals in 18 populations were collected (excluding population PSET). Linear mixedeffects models were conducted on the four variables as described in the "Measurements and Analyses of Leaf Morphology" section.

| Macrohabitat differentiation
The macrohabitat axes of bioclimatic data were retrieved based on range-wide occurrences at the coarse scale. In addition to our field sampling localities (19 records,

| Geographical cline analysis
To explore whether endogenous or exogenous selection is primarily responsible for maintaining the hybrid zone, one-dimensional (1D) cline analysis was performed based on both genetic data (frequency of O. leptophylla cpDNA haplotype and microsatellite Q K2 scores) and environmental data (log-transformed slope values and macroclimatic PC2 scores). Morphological data (leaf-shape frequency and leaf morphology PC2 scores) were also used to fit geographical cline to test for any effect of divergent selection or morphological introgression. Three-model comparison were conducted to fit either allele frequency data or quantitative trait data (for details of these models see Brumfield et al., 2001), and the best-fit model with the lowest Akaike information criterion (AIC) value was selected to estimate cline center and width.
Third, the program Cfit7 (Gay et al., 2008) was used to test for cline coincidence and concordance. This program allows two kinds of models compared: one in which different clines constrained to a common center and/or width, and one in which the cline shapes were allowed to vary independently. The constrained model was treated as significantly better than the others using AIC (i.e., the convention of selecting one that has an AIC value smaller by more than two AIC points).

| Population genetics
Among the 19 locations sampled, 33 out of 1045 (3.2%) withinpopulation genotypic disequilibria were detected prior to correction for multiple testing, but were not specific for any pairwise loci or population. After Benjamini-Hochberg correction, only one pair of linkage disequilibrium in population PSEV remained significant (Table S6). Significant heterozygote deficits were observed at 1 to 10 loci in each population, and in total 138 out of 209 locus-specific tests was significant (

| Parent and hybrid assignment
According to the STRUCTURE analysis, the optimal K value was 2 ( Figure S2), with the two parental species formed distinct genetic groups and the hybrid zone showed genetic admixture (Figure 3).  (Table 3)  TA B L E 2 Descriptive statistics of genetic variability based on 11 nuclear microsatellite loci and five cpDNA intergenic regions.

F I G U R E 3 STRUCTURE results for 19
Oxytropis populations based on microsatellite dataset (K = 2-3). The small black lines separate populations, named at the bottom of the graph. Each individual is a small bar with color coded according to probability of clustering in a particular group. Populations are ordered west to east for each species. greater than a 9.5% chance of being classified as F 1 s ( Table S4). This analysis struggled to assign individuals to specific classes with any confidence, which is probably due to the lack of species-diagnostic loci (Anderson & Thompson, 2002), and could also be an expected outcome when hybrid derivatives of the third generation or beyond, and other complex intermediate hybrid derivatives are present (Yan et al., 2017).

| Chloroplast DNA sequence analyses
The concatenated cpDNA dataset had a total alignment length of 3918 bp. Across all 261 sequences, 57 haplotypes were identified, characterized by 80 variable sites, of which 45 were parsimony informative.
The 197 individuals from putatively "pure" populations of two parental species were fixed for distinct haplotypes (Table 3), with H1-H51 in O. diversifolia and H54-H57 in O. leptophylla, respectively ( Figure S3). Two haplotypes, H52 and H53, were only found in the hybrid zone, but they were derivatives of O. diversifolia haplotypes H42 and H6, respectively ( Figure S3). Consequently, they can be O. diversifolia haplotypes miss detected in putatively "pure" populations. Significant cytonuclear disequilibrium was detected in both population G and HUA (p < .01, Table S6). The genetic diversity estimates of S and π were highest in TA B L E 3 Number of parental and hybrid individuals assigned based on nuclear microsatellite data  Figure S3). Of these, 29 individuals were classified as hybrids with nuclear analyses (Table S4) ratio, not differing significantly from the 1:1 expectation for no gender bias (χ 2 = 0.86, df = 1, p = .35; Table 3).

| Microhabitat differentiation
We found that values of slope (°) in the hybrid zone (10.1 ± 6.41, mean ± SD) were intermediate between the two parental species (4.95 ± 3.63 in O. diversifolia and 12.9 ± 9.74 in O. leptophylla, respectively) (p < .01; Figure 5a, Table S9). However, the percent of total vegetation cover, the percent of rocky ground, and the percent of bare ground did not differ among categories (p = .15-.69; Table S9).

| Macrohabitat differentiation
Seven out of the 19 bioclimatic variables showed significant differentiation, with localities of the hybrid zone possessed intermediate values between those of the two parental species (Table S10) luvic Kastanozems (KSl) ( Figure S5b). Similarly, the associated plant species in 19 localities that we surveyed were coupled to the gradual change of both macroclimatic and soil characteristics (Table S11).

| Geographical cline analysis
We plotted the maximum-likelihood cline for the best-fit model for each morphological, genetic, and environmental variable (Figure 6a-f). All variables, except for macroclimatic PC2 score, ex-

| DISCUSS ION
Hybrid zones generated through secondary contact between closely related species have been widely highlighted for the interest in understanding the evolutionary processes (Barton & Hewitt, 1985;Harrison, 1993). However, the number of plant hybrid zones detected and studied genetically was rather limited (reviewed in Abbott, 2017).
In the present study, we integrated morphological, molecular, and ecological analyses to investigate the hybrid zone generated between

| Endogenous selection primarily responsible for maintaining the narrow hybrid zone
Hybrid zones are often formed after the secondary contact of two divergent parental forms, and some previous definitions restricted the term hybrid zone to narrow clines (Barton & Hewitt, 1985), although mosaic hybrid zones can also be common (Abbott, 2017).
In this study, the two Oxytropis species are separated by a narrow hybrid zone that can be maintained by a balance between dispersal and divergent selection. This inference is based on the shape of the phenotypic and genetic transition (i.e., sigmoidal clines), the narrow width relative to the dispersal ability (Wang et al., 2021), and the elevated variance at the center (Barton & Hewitt, 1985;Stankowski et al., 2021). However, solely from the sigmoidal cline, it is difficult to distinguish between selection against hybrids (i.e., "endogenous" selection) and adaptation to different environments (i.e., "exogenous" selection), because endogenous and exogenous selection can generate clines of similar shape (Barton & Gale, 1993;Barton & Hewitt, 1989;Kruuk et al., 1999).
Taken our results together, we suggest that endogenous selection can be more important. First, although the genetic clines coincide with the environmental cline, the former is significantly narrower than the latter. This is particularly evident for the comparison between genetic clines and the macroclimatic gradient. The observed pattern is consistent with the prediction of tension zone model, in which endogenous selection can create a sufficiently strong barrier to gene flow independent of habitat (Barton & Hewitt, 1985, 1989. Second, we investigated clinal variation across the hybrid zone in a variety of characters, either morphological or genetical, which show close concordance as parallel gradients. If each allele was to respond directly to the environment, such geographic concordance would not be expected (Barton & Hewitt, 1985). Finally, the evidence for endogenous selection against hybrids is usually in the form of hybrid inviability, sterility or abnormality in the field or laboratory (Stankowski et al., 2021). Due to the perennial feature and field work limitation, we are not able to provide a direct measure of hybrid inviability or sterility (but see Abbott, 2008 andYan et al., 2019 for example). However, we did observe abnormality of leaf characteristics in the hybrid zone (Table S4).
Because tension zones can be maintained independent of environmental conditions, they can move from place to place. There are two main explanations for the position of hybrid zone. The first hypothesis is the "low-density troughs". Environmental features like valleys, rivers, ridges and wood cause density to vary over orders of magnitude, and moving hybrid zones have the interesting property of tending to rest in regions of low density (Hewitt, 1988). The alternative hypothesis suggests that tension zones have a tendency to become trapped by, and therefore to coincide with, exogenous barriers due to ecological selection (Bierne et al., 2011). In our study, none of the clines are positioned in an area of low density. Instead, genetic clines are coupled to ecological transitions. Theoretical simulation indicates that coupling is more likely with moderate exogenous selection, as it results in very wide cline (e.g., latitudinal clines) able to trap the endogenous tension zone from a remote distance (Bierne et al., 2011). In our case, the microhabitat cline is wider than genetic clines, and the macroclimatic transition is even wider along the longitude. Both of them can produce moderate exogenous selection to trap the tension zone.

| Potential effects of hybridization
A pronounced feature of the hybrid zone detected in our study is the bimodal pattern of genotypic distribution. Historical literature review suggests that both bimodal and unimodal hybrid zones have been shown to exhibit postzygotic isolating mechanisms, which may be either intrinsic or extrinsic. However, only bimodal zones are associated strongly with prezygotic barriers (Jiggins & Mallet, 2000). acterization of these hybrids using species-specific nuclear markers should be warranted (e.g., SNPs).
Our data also indicated potential cytoplasmic introgression.
We found that two rare cpDNA haplotypes of O. leptophylla, i.e., H56 and H57 detected in the allopatric population ZHONG, are distinct from major haplotypes of this species but close to O. diversifolia haplotypes ( Figure S3), while the two individuals possessing those haplotypes are still of O. leptophylla morphology and nuclear germplasm. This may be due to an ancient "chloroplast capture" event through hybridization and introgression, which frequently happens in plants (Rieseberg, 1995). Chloroplast capture is probably one of the most common causes of phylogenetic incongruence between nuclear and plastid gene tree (Wendel & Doyle, 1998).
It has been reported in phylogenetic studies involving deep relationships among different genera (e.g., Liu et al., 2017), and also detected at lower taxonomic ranks in population and/or specieslevel studies (e.g., Wendel et al., 1991). Moreover, in the hybrid zone we found six individuals with the nuclear germplasm of O.
diversifolia but cpDNA haplotypes of O. leptophylla and vice versa ( Table S4), suggesting that there might be ongoing cytoplasmic introgression.
Introgression is emerging as an important source of novel genetic variation, alongside with standing variation and mutation, and it can be adaptive when such introgressed alleles are maintained by natural selection (termed "adaptive introgression", Abbott et al., 2013;Suarez-Gonzalez et al., 2018). A previous study revealed that the clinal variation in leaf shape in O. diversifolia can be attribute to natural selection imposed by spatially varying macroclimatic factors (Wang et al., 2021). It would be interesting to test whether leaf-shape polymorphism in O. diversifolia was obtained through hybridization and introgression from O. leptophylla to the 1-leaflet phenotype. The detection of nuclear introgression requires makers that are strongly species-specific, and again, microsatellites used in our study is not satisfying. Several studies have applied next-generation sequencing and revealed that introgression transfers adaptive alleles. For example, Arabidopsis arenosa became adapted to drought and toxic levels of soil minerals in serpentine soil through introgression of alleles from A. lyrata (Arnold et al., 2016). Similarly, introgression from Cupressus gigantea into a related species C. duclouxiana likely aided the latter species to extend its range by colonizing cooler and drier mountain habitats during postglacial periods (Ma et al., 2019).

| CON CLUS ION
In summary, we have found evidence of bimodal hybrid zone connecting two parapatric Oxytropis species; the geographic cline analysis indicated that the hybrid zone can be maintained by a balance between dispersal and selection, with the latter primarily being endogenous. We further conclude that the probable outcome of hybridization is introgression. A potential hypothesis to be tested is adaptive introgression, i.e., whether leaf-shape polymorphism in O.
diversifolia was obtained through introgression from O. leptophylla.

ACK N OWLED G M ENTS
We thank Liang Xiao, Chun Su, Yuan Lu, and Lu-Lu Xun for their help with the field sampling, and Han Yang for her assistance with lab work. We are grateful to Professor Mitch Cruzan and three anony-

CO N FLI C T O F I NTE R E S T
The authors declare no conflict of interest.